www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/APfilter.m

    function [ b ] = APfilter( a,w,SCmask,SCit )
%APFILTER Adaptive Phase Filter
%  Syntax:
%  [ b ] = APfilter( a,w,SCmask,SCit )
%  
%  a is a phase map
%  w is the size of the filter window
%  SCmask is mask for Sin/Cos filter
%  SCit is the iteration time for Sin/Cos filter
%  ------------------------------------------------------------------------
%  ------------------------------------------------------------------------
%  Reference:
%  [1]
%  Francisco Palacios, Edison Goncalves, Jorge Recardo, etal. Adaptive
%  filter to improve the performance of phase-unwrapping in digital 
%  holography. Optics Communications. 238,2004:245–251
%  ------------------------------------------------------------------------
a=SCfilter( SCmask,a,SCit );
[M,N]=size(a);
Gamma=nint((a(1:M-1,2:N)-a(1:M-1,1:N-1))/2/pi);       % d=A
Gamma=Gamma+nint((a(2:M,2:N)-a(1:M-1,2:N))/2/pi);     % d=A+B
Gamma=Gamma+nint((a(2:M,1:N-1)-a(2:M,2:N))/2/pi);     % d=A+B+C
Gamma=Gamma+nint((a(1:M-1,1:N-1)-a(2:M,1:N-1))/2/pi); % d=A+B+C+D
Gamma=[Gamma,zeros(M-1,1)];
Gamma=[Gamma;zeros(1,N)];
T=conv2(Gamma,ones(w(1),w(2))/w(1)/w(2),'same');  % NDP
clear Gamma
%
if rem(w(1),2)==0
    om=w(1)/2+1;
    c=[zeros(w(1)/2,N);a;zeros(w(1)/2-1,N)];
else
    om=(w(1)+1)/2;
    c=[zeros((w(1)-1)/2,N);a;zeros((w(1)-1)/2,N)];
end
if rem(w(2),2)==0
    on=w(2)/2+1;
    c=[zeros(M+w(1)-1,w(2)/2),c,zeros(M+w(1)-1,w(2)/2-1)];
else
    on=(w(2)+1)/2;
    c=[zeros(M+w(1)-1,(w(2)-1)/2),c,zeros(M+w(1)-1,(w(2)-1)/2)];
end
[X,Y]=meshgrid(1:w(2),1:w(1));
preg=((X-on).^2+(Y-om).^2)/2;
d=zeros(M,N);
for m=1:M
    for n=1:N
        g=exp(preg/(T(m,n)+1).^2);
        g=g/sum(sum(g));
        %d(m,n)=sum(sum(c(m:m+w(1)-1,n:n+w(2)-1).*g));
        d(m,n)=angle(sum(sum(exp(i*c(m:m+w(1)-1,n:n+w(2)-1)).*g)));
    end
end
%
e=zeros(M+w(1)-1,N+w(2)-1);
e(om:om+M-1,on:on+N-1)=d;
clear d
b=zeros(M,N);
for m=1:M
    for n=1:N
        g=exp(preg/(T(m,n)+1).^2);
        f=T(m,n)*abs(e(m:m+w(1)-1,n:n+w(2)-1)-e(m+om-1,n+on-1))+1;
        h=g.*f;
        h=h/sum(sum(h));
        %b(m,n)=sum(sum(c(m:m+w(1)-1,n:n+w(2)-1).*h));
        b(m,n)=angle(sum(sum(exp(i*c(m:m+w(1)-1,n:n+w(2)-1)).*h)));
    end
end

function b=nint(a)
%  round numbers between -1 and 1 to nearest integer
%  0.5 and -0.5 are rounded to 0
b=sign(a)-round(sign(a)-a);